### Do a quick plot of density of accuCLIM05 vs BIOCLIM05

### out directory

out.dir = '/home1/99/jc152199/Paluma/'

### load library

library('SDMTools')

## Read in ASCII's

ac = read.asc.gz('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/microCLIM/microCLIM_05.asc.gz')
fivek = read.asc('/home1/99/jc152199/MicroclimateStatisticalDownscale/5kASCII/AnalysisMask5kResolution/ASCII/projected.asc')

### Remove NA's to determine density

ac2d = ac[which(is.finite(ac)==T)]

### Determine the 5k points from accCLIM

pos = as.data.frame(which(is.finite(fivek), arr.ind = T))

### Get lat/longs

pos$lat = getXYcoords(fivek)$x[pos$row]
pos$long = getXYcoords(fivek)$y[pos$col]

### Intersect with ac

ac5k2d = extract.data(cbind(pos$lat,pos$long),ac)

### Remove NA's

ac5k2d = ac5k2d[which(is.finite(ac5k2d)==T)]

# Determine density

ad250 = density(ac2d)
ad5k = density(ac5k2d)

### Xlims

xlims = range(c(ad250$x))
xlims = round(c(xlims[1]-1,xlims[2]+1),1)


### Y lims

ylims = c(0,max(c(ad250$y, ad5k$y)))
ylims[2] = round(ylims[2]+.01,2)

### Open .png device driver

png(paste(out.dir,'05overlay.png',sep=''),height=18,width=18, res=500, units='cm')


### Configure plot space

plot(ad250$x,ad250$y, main='Mean Max Temp', xlab='MMTofWP', ylab='Density', type='n', xlim=xlims, ylim=ylims)

### Plot lines

points(ad250$x,ad250$y, col='blue', type='l')
points(ad5k$x,ad5k$y, col='red', type='l')

### Add a legend

legend('topleft',legend = c('250 m','5 k'), text.col = c('blue','red'), bty='n')

### Close device driver

dev.off()

#### Next idea #####################################################################################

png(paste(out.dir,'hist3.png',sep=''), units='cm',width = 18, height=9, res=1000)
par(mfrow=c(1,2))
hist(ac2d, xlim = range(ac2d))
hist(ac5k2d, xlim = range(ac2d))
dev.off()

png(paste(out.dir,'hist2.png',sep=''))
hist(ac5k2d)
dev.off()

